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Abstract 

The probabihty distribution p{l) of an atom to return to a step at distance I from the de- 
tachment site, with a random walk in between, is exactly enumerated. In particular, we 
study the dependence of p{l) on step roughness, presence of other reflecting or adsorbing 
steps, interaction between steps and the diffusing atom, as well as concentration of defects 
on the terrace neighbouring the step. Applying Monte Carlo techniques, the time evolution 
of equilibrium step fluctuations is computed for specific forms of return probabilities. Re- 
sults are compared to previous theoretical and experimental findings. 
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1. Introduction 



In recent years, the dynamics of steps on crystal surfaces has attracted much interest, 
both experimentally [1-6] and theoretically [7-11]. Experimentally, equilibrium step fluctu- 
ations of isolated steps as well as steps on vicinal surfaces have been studied extensively, 
following pioneering STM measurements of Au(llO) [1] and Cu(lln) [2] surfaces. The step 
fluctuations are quantified by the correlation function G{t)— {{h{i,t) — /i(i,0))^), where 
h{i, t) denotes the position (or displacement) h of the step at site i and time t. Typically, 
the experimental data are described by a power law, G{t) oc t^, with the dynamic exponent 
7 being in between about 1/4 and about 1/2. Theoretically, distinct atomic mechanisms 
driving the step dynamics have been identified, leading, indeed, at long times, to power 
laws, with 7 being 1/4, 1/3, and 1/2 in the limiting cases of step diffusion, terrace diffusion, 
and evaporation-condensation kinetics, respectively. These scenarios have been found in 
Langevin descriptions [7,9,10] and confirmed, partly, in simulations on SOS models [7,11]. 

In this article, step fiuctuations will be analysed in the framework of random walks of atoms 
detaching from the (reference) step, diffusing on the neighbouring terrace, and attaching 
again at the step at distance / (in units of the lattice spacing on the surface) from the 
detachment site. The corresponding return probability distribution p{l) is determined for 
various situations, to analyse the influence of the step roughness, of the presence of another 
reflecting or adsorbing step separated by d lattice spacings from the reference step, of the 
interaction of the diffusing atom with the steps due to, e.g., elastic forces, as well as of the 
concentration of reflecting or adsorbing defects on the terrace. The case of a straight step 
with no obstacles on the terrace and without interactions between the atom and the step 
has been considered before, solving the related continuum diffusion equation, in the context 
of a Langevin theory, where p{l) is called diffusion kernel [9]. Here, results on p{l) will be 
then incorporated in Monte Carlo simulations of the step fluctuations. We shall present 
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simulational results for special choices of return probabilities connecting detachment and 
attachment step sites, including the three limiting cases mentioned above. One of the main 
points will be to elucidate the relation between those probability distributions and the step 
fluctuations characterised by G{t). Indeed, deviations from the simple power law for G{t) 
are found, which may play a relevant role in interpreting experimental observations. 

The article is organized as follows. In the next section, the model used to describe step 
fluctuations is introduced, and the methods applied in the analysis, exact enumeration of 
random walks and Monte Carlo simulations, are introduced and outlined. Results on the 
return probability distribution p{l) are given in section 3, followed by a discussion of simu- 
lational findings on G{t). A brief summary will conclude the paper. 

2. Model and methods 

Let us consider a (reference) surface step of monoatomic height with 2L + 1 sites, 
i — — L, —L + 1, 0, ...L — 1, L. The positions of the step atoms, at time t, are denoted 
by h{i,t). Setting the lattice spacing equal to one, h is an integer. The step is perfectly 
straight, if h is constant. The neighbouring terrace, consisting of sites (i, j) on a square 
lattice, may be bordered by a straight step opposite to the reference step, separated by d 
lattice spacings, see Fig. 1. The opposite step, with 2L-I- 1 sites, is supposed to either reflect 
(variant s,.) or absorb (s^) an adatom diffusing on the terrace, mimicing either a descending 
step with a high Schwoebel-Ehrlich barrier or a rising step with a large sticking coefficient. 

Specific boundary conditions may be used at the ends of the reference step, i = —L and 
i — L. For instance, the ends may be pinned or subject to periodic boundary conditions. 
Alternatively, the boundaries of the terrace perpendicular to the reference and opposite steps 
may be, e.g., chosen as straight reflecting (variant tr) or absorbing (ta) straight steps with 
d sites. 
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Step fluctuations result from the detachment of an atom from the reference step at site 
i, h — > h — 1, and attachment at step site i + I, h — > h + 1. The rate of this process 
depends on the energies involved and temperature. One may distinguish three limiting cases, 
(a) / = ±1, corresponding to step diffusion; (b) uncorrelated detachment and attachment 
sites, corresponding to evaporation-condensation kinetics; and (c) initial and final sites 
being connected by a random walk of the diffusing atom on the terrace neighbouring the 
step, corresponding to terrace diffusion. Several types of terrace diffusion are possible, 
reflecting various constraints on the random walk, see below. The time dependence of the 
step fluctuations may be quantifled by the correlation function 

G{t) = {{h{t,t)-h{t,Q)y) (1) 

The brackets denote a thermal average. 

We applied two methods in our study: (i) exact enumeration of random walks [12], to inves- 
tigate different kinds of terrace diffusion, and (ii) Monte Carlo simulations [13], to compute 
G{t). 

(i) In enumerating random walks, we calculate the return probability distribution p{l) of 
an atom detaching from site io of the reference step onto the terrace (the step positions 
being then flxed once and for all, /io(^) — h{i,t)) and attaching at step site iq + 1. In other 
words, the atom starts its random walk at terrace site (io, /lo(^o) +2), and one determines the 
probability that it reaches the site (io + I, ho{io + /) + 1) first among all possible absorption 
sites at the reference step. The terrace is supposed to be bordered by either reflecting or 
absorbing steps perpendicular and opposite to the reference step. Obviously, one has 

Ep(0 < 1 (2) 

I 

where the sum runs over all step sites. The identity holds when the boundary steps are 
reflecting. Usually, the starting point of the random walk is chosen at the center of the step, 
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The exact enumeration of p{l) is based on calculating the probability Prw(^,J, ^) of the 
diffusing atom to be at terrace site (i, j) in the n-th move of the random walk [12,14]. The 
total probability, summing PTw{hj,n) over all sites, is conserved and equal to one. p^w is 
related to the probability to jump from the neighbouring site (i', j') to site (i, j), w{i', j', i, j), 

by 

Prw(«, j, n)^Yl ^(^'' ^' j)Pr^i^'^ /' n-l)+ Ws(i, j)pr^(i, j, n - 1) (3) 
i'j' 

where Ws{i,j) is the probability to stay at terrace site For an absorption site, at the 

reference or a boundary step, one has Ws{ij) — 1 and w{i, , j')—0, while for the other 
terrace sites Ws{i,j) = 0. 

On finite terraces, pj^{i,j,n) can be calculated in a straightforward fashion for a variety 
of jump probabilities and boundary conditions, with the bookkeeping done by a computer 
programme (exact analytic results are known only for rather few special cases [15]). The 
return probability distribution p{l) is given by 

p{l) ^Prwiio + l,hoiio + l) + 'i-,n) (4) 

considering indefinitely long random walks, n — > oo. In practice, the length of the walk, n, 
depends on the convergence rate of p^w- We studied steps of length L < 600, with the width 
of the terrace d < 600. Typically, the random walks ended when the total probability of 
finding the atom on any non- absorption terrace site was smaller than 10~^. 

(ii) In our Monte Carlo simulations an atom at step site i is moved to site i + l with a (return) 
probability Pda(0- The form of Pda(0 is motivated by the findings on the random walks, p{l). 
The move is then accepted, as usual, with a rate determined by the Boltzmann factor of the 
associated energy change SE, exp{—SE/kBT), where ks is the Boltzmann constant and T 



the temperature [13]. The energies at the step are assumed to be given by the number and 
depth of the kinks [16], as it is the case in the standard SOS model [17] with the Hamiltonian 

n = cj2\h{i,t)-h{j,t)\ (5) 

where the sum runs over neighbouring step sites j = i ± 1.- To speed up the simula- 
tions, we used an algorithm with a dynamic time assignment [13,18]. 

The time t, elapsed during the simulation, is measured in terms of Monte Carlo attempts 
per step site pair (MCA), i.e. in one time unit one has tried to interchange, on average, one 
atom between each two step sites (the Monte Carlo time scale is linearly related to the real 
time) . This interpretation corresponds to the situation where the time spent by the diffusing 
atom on the terrace is negligibly small compared to the mean time spent at the step. More 
realistic approaches using, e.g., kinetic Monte Carlo simulations for the surface dynamics, 
may suffer from other shortcomings, such as an ambiguity in the step position. Note that 
the average position of the entire step, J2h{i,t)/{2L + 1), does obviously not depend on 

i 

time. 

To compute the step fluctuations G{t), see eq. (1), one has to define a starting time, t = 0. 
Various choices are possible, corresponding to various initial step configurations, including 
straight and equilibrated steps. Perhaps closest to experiments, one may average over an en- 
semble of thermalized step configurations, possibly during a single Monte Carlo run (which 
will be called, in the following, dynamic averaging). In addition, one may impose different 
boundary conditions, for instance, by pinning the end positions of the step, at i = 1 and 
2L + 1, or applying periodic boundary conditions. In the later case, boundary effects are 
less severe. 

Typically, we studied steps of length 2L + 1 < 256, with periodic boundary conditions. 
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3. Random walks 



Using the exact enumeration approach, the return probabihty distribution p{l) has been 
calculated for four distinct cases, mimicing possible constraints on the terrace diffusion: (a) 
non-perturbed random walk, w{i' , j' ,i, j) = 1/4, with a straight reference step, see Fig. 2; 
(b) diffusing atom in an external potential, V = A/y"^, y being the distance between the 
(straight) steps and the atom, due to, e.g., elastic interactions [19], see Fig.3. Note that due 
to the potential V the jump probability w oc exp(— 51//A;bT), with 5V = V{i',j') — V{i,j), 
becomes anisotropic, favouring hops away from the step; (c) unbiased diffusion, w — 1/4, 
on a terrace with quenched absorbing single-site defects of concentration c, with a straight 
reference step, see Fig.4; and (d) non-perturbed random walk with a rough reference step 
(the roughness mimics a thermalized reference step). 

In each case, the wandering atom started at the center of the reference step, at terrace site 
(0, /io(0) + 2). The three boundary steps were assumed to be straight and either reflecting 
or absorbing. For instance, 's^ '■ ta refers to a situation with a reflecting opposite step and 
absorbing terrace boundary steps perpendicular to it. 

In general, one may distinguish four regimes, in which p{l) exhibits different characteristic 
properties, depending on the return distance I. At very short distances of a few lattice 
spacings, I < Iq, p{l) falls off rapidly, the concrete form being determined by details of the 
perturbations (step roughness, concentration of defects and strength A of the elastic inter- 
action) . 

At Iq <^ I <^ d, p{l) acquires typically a power-law behaviour, p{l) oc Z". To analyse that 
regime, one may calculate the effective exponent 

a^s{l)^d\np{l)/d\nl (6) 
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being constant, a^s — if, indeed, the deacy of p{l) follows a simple power-law. Increasing 
the return distance I furthermore, the effect of the boundaries on the random walk shows 
up. Assuming L ^ d, first the presence of the opposite (vicinal) step affects the return 
probability, leading to a characteristic exponential decay of p{l), p oc exp(— aZ), both for 
reflecting, Sr, and absorbing, Sa, steps. Finally, the distribution p{l) will be modifled due to 
the perpendicular terrace boundaries, tr or ta- 

As illustrated in Fig. 2 for the non-perturbed case, cteff tends to approach a — —2, in 
agreement with the solution of the corresponding continuum diffusion equation [9]. Con- 
sequently, one expects p{l) — pol~^ at sufficiently large return distances I for indeflnitely 
long isolated straight steps (note that a is expected to determine the value of the dynamic 
exponent 7 describing the time dependence of the step fluctuations G{t) [9], see below). 
Actually, neither the interaction of the diffusing atom with the steps nor the step roughness 
seem to affect that value of a. In fact, for a rough reference step, aes{l) follows closely 
the form for a straight step, / > Iq [14]. Applying the potential V = A/y"^, we observe a 
systematic dependence of the proportionality factor po on A, as shown in Fig. 3. po is found 
to increase exponentially with the interaction strength A, pq oc exp{r]A), with rj ~ 1.23, at 
least for the strengths we considered, < A < 1. As depicted in Fig. 4, our findings do 
not rule out that the value of a may, however, depend on the concentration c of defects on 
the terrace. Over an appreciable range of return distances I one may notice a plateau-like 
behaviour in a^s at a value slightly, but definitely smaller than -2. E.g., at c = 0.2, the 
plateau-like behaviour occurs at about —2.2; at c = 0.1, the effective exponent ccefr seems 
to settle for some distances I at a value below, but closer to -2. However, the data do 
not allow to rule out the possibility that a = —2 may be approached, now from below, at 
larger return distances I. Of course, it would be desirable to consider longer steps and larger 
terraces to study further this aspect. Such computations would be extremely computer time 
consuming, because one has to average over a large number of defect concentrations to get 
meaningful data (typically, we averaged over about 10 000 realizations). 
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With increasing /, the return probabihty distribution p{l) of the diffusing atom approaches 
an exponential decay, p{l) oc exp{—gl) in all four cases, {a)-{d). Accordingly, the expo- 
nential form describes p{l) for large return distances I in the presence of a neighbouring, 
descending or rising, step, separated by d lattice spacings. For the unbiased random walk 
with straight steps the coefficient g, as follows from exact enumeration, agrees, to a high 
degree of accuracy, with the one obtained from the continuum diffusion equation [9] , namely 
g — 7r/rf for an absorbing opposite step, and g = Tr/2d ior a. reflecting, vicinal step. The 
various perturbations seem to have little effect: The decay form of remains exponential, 
and even the coefficient g seems to be rather robust.- In the unbiased case, (a), we found 
that p{l, d) seems to scale, at d ^ 1, as p{l, d) = Ps{l/d)/d'^, both in the power-law and the 
exponential regime [14]. 

Finally, with / getting closer to L, the perpendicular terrace boundaries become relevant. As 
expected, p{l) is lower when the boundaries are absorbing than in the reflecting situation. 
However, we did not explore this region in much detail, because those boundary conditions 
are somewhat artiflcial and extremely long random walks, see eq. (4), would be needed for 
quantitatively correct results. 

4. Monte Carlo simulations 

We simulated the fluctuations G{t) of steps of length M = 2L + 1, with their ends being 
connected by periodic boundary conditions. An atom is detached at the randomly chosen 
site i and attached at site i + /, see Fig. 1. The (return) probability of selecting site i + / 
is denoted by Pda(0- The probability of accepting such an elementary move is given by the 
Boltzmann factor of the energy change, see eq. (5), associated with this process. 

We studied mainly flve distinct cases, with the choice of pda(0 being motivated by results on 
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the return probability distribution p{l) of a random walk and by general physical considera- 
tions. In particular, we considered Pda(0 being (i) constant, corresponding to evaporation- 
condensation kinetics. Obviously, the mean step position is conserved in each move. We 
also simulated evaporation-condensation dynamics by choosing between detachment and 
attachment completely randomly, thus conserving the mean step position only in the time 
average; (ii) 1/2 for / = ±1 and elsewhere, describing step diffusion; (iii) proportional 
to 1/P (the prefactor is determined by requiring that Pda(0 normalized: Z^Pda(0 — 1)> 
as one may expect for fluctuations of an isolated step due to terrace diffusion, neglecting 
deviations from that form at short distances / between attachment and detachment sites; 
(iv) proportional to exp(— Z), monitoring the possible effect of a neighbouring step on G{t); 
and (v) proportional to to study the sensitivity of G{t) on the value of a in Pda(0 °^ 

For indefinitely long steps, one may argue that G{t) increases with time in the form of 
a power-law G{t) oc t"^ [7-10], with the dynamic exponent 7 depending on the atomistic 
mechanism governing the fluctuations. To determine 7 in the various cases we simulated, 
one may calculate the effective dynamic exponent 

-f^s(t) ^ din G{t)/ dint (7) 

with 7efr(^) approaching 7 at large times and long steps. 

Results on 7efr(t) for each of the five situations are depicted in Fig. 5, at fixed step length, 
M — 256, and fixed temperature ksT/e — 1. We used dynamic averaging over an ensem- 
ble of successive, equilibrated initial configurations. Evaporation-condensation kinetics was 
simulated by random attachment and detachment processes, i.e. without conserving the 
average step position in each Monte Carlo move. 

At very early times, the step fiuctuations are always diffusive, with G{t) = Cdt, corresponding 
to 7eff = 1. The diffusion coefficient q depends rather weakly on the transport mechanism, 
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with the exception of the evaporation-condensation dynamics having a significantly larger 
diffusion coefficient. 

Due to the rigidity of the step, the step meandering then slows down, with 7eff approaching 
a plateau located at about 1/2, 1/3, or 1/4, depending on the form of the return probability, 
Pda, see Fig. 5. Actually, the plateau at 7efr ~ 1/2, characterising evaporation-condensation, 
is reached most quickly. The plateau at 1/3 is reahsed for terrace diffusion with an isolated 
step, Pda(0 ^ 1/^^- The largest time is needed to approach the plateau at 1/4, character- 
ising step diffusion as well as terrace diffusion with rapidly decaying return probabilities, 
Pda(0 oc 1/^^ and oc exp(— Z). The last situation corresponds to terrace diffusion with pairs 
of steps. One may argue that the plateau signals the asymptotic power law increase of 
G{t) oc P for large times and indefinitely long steps. Indeed, the plateau values have been 
obtained before in Langevin descriptions [7-9] for infinitely extended steps (for —3 < a < —2, 
one may expect a continuously varying value, 7 = l/(— a + 1) [9]). 

Combining our results on the random walks and simulations, we conclude that the exponent 
7 = 1/3 for an isolated step is robust against various perturbations of terrace diffusion, with 
the possible exception of defects on the terrace. 

The estimates of the plateau values of 7eff are confirmed by a Fourier analysis of the step 
fluctuations h{i,t). They can be written in the form 

M/2 

h{i, t) = Y^ (afc(t) sin(27rA;i/M) + bk cos(27rA;i/M)) (8) 
fc=i 

The fluctuation modes may be described by 

Gk{t) = {{ak{t) - akiO) f + (hit) - hiO) f) (9) 

taking into account phase shifts. Based on Langevin descriptions, Gfc(t) is expected to 
converge rapidly towards equilibrium Gk{oo) [7-9], 
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Gk{t) = Gk{oo){l - exp(-7fei)) (10) 

in agreement with the simulational data. More specificly, for small wavenumbers, Ik follows 
closely the form Ik oc k^, with x = I/7, in accordance with the Langevin theory However, Ik 
shows pronounced deviations from the power-law behaviour at larger values of /c, as shown 
in Fig. 6. This feature is not described by the Langevin theory [7]. 

At later times, after having passed the plateau, the effective dynamic exponent, 7efr, will 
eventually go to zero (as may be easily seen for shorter steps [11]) or, in the case of 
evaporation-condensation, to one, due to the constraint of conserving, or not, the aver- 
age step position, see Fig. 5. For evaporation-condensation, the related diffusion coefficient 
is rather small and depends on the step length. In any event, the late-time behaviour reflects 
the finite length of the step. 

The time dependence of %s, as depicted in Fig. 5, implies that an average effective ex- 
ponent, 7a, as obtained from a log-log plot of G{t) in a given time range, may vary with 
experimental parameters when one is not in the truely asymptotic regime. In particular, the 
plateau in ^y^s may not yet have been reached or finite size effects may already matter. Note 
that, for instance, the crossover from the diffusive short-time behaviour to the subdiffusive 
step motion at later times is temperature dependent, thereby causing possibly a tempera- 
ture dependence in 7a. Similarly, the presence of a neighbouring step is expected to affect 
7a. In general, much care is needed in identifying and disentangling the numerous possible 
crossover effects. 

5. Summary 

We studied equilibrium step fluctuations in a somewhat idealized way by enumerating 
random walks on terraces with various constraints and using Monte Carlo techniques. 
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In the simulations, different atomistic mecfianisms determining tfie equilibrium step fluctu- 
ations may be mimiced by special choices of the return probability for atoms attaching at 
distance I from the detachment site of the step. The form of the return probability at large 
distances / determines the time dependence of the step fluctuations G{t) at late times t for 
long steps. In particular, we simulated evaporation-condensation kinetics, step diffusion and 
(perfect) terrace diffusion. For isolated steps, G{t) is confirmed to increase with a power 
law G{t) oc n, with 7 =1/2, 1/3, and 1/4, respectively. In the case of terrace diffusion for 
pairs of steps, the dynamic exponent 7 approaches the value of step diffusion, 7 = 1/4. The 
simulational observations agree with and refine predictions of Langevin theory.- Similarly, 
a Fourier analysis of the step fluctuations, driven by the different atomic mechanisms, con- 
firms and refines previous Langevin descriptions. 

In general, crossover phenomena may mask the asymptotic behaviour of G{t). Such phenom- 
ena can be caused by several reasons, including change from diffusive to subdiffusive step 
motions at early times, effect of the finite step length, defects on the terrace, and influence 
of neighbouring steps. As a result, the average effective dynamic exponent 7^, as usually 
obtained from measurements, may vary with experimental parameters like temperature. 

From the exact enumeration of random walks, the decay of the return probability distribu- 
tion p{l) oc for unbiased terrace diffusion with straight steps is found to be robust against 
step roughness and (elastic) interactions between the diffusing atom and the neighbouring 
step. This fact implies the robustness of the value 7 = 1/3 for isolated steps. Moreover, for 
pairs of steps, the exponential decay of p{l) persists in the presence of those perturbations, 
implying the robustness of 7 = 1/4 in that case. 
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Figure Captions 

Fig. 1: Geometry of the surface, showing the reference step and a descending (s^) or rising 
(sa) opposite step. 

Fig. 2: Effective exponent aes of the return probabihty distribution p{l) for the unbiased 
random walk with a straight reference step (L = 600, d = 600). Absorbing terrace 
boundaries, ta, were used. 

Fig. 3: Return probability distribution p{l) for random walks in a repulsive external potential 
V — Ay~^ at ksT — 1, ior L — 600, d — 600, and with : ta boundary conditions. 

Fig. 4: Effective exponent of p{l) for random walks with defects on the terrace {d — 50, 
L = 100, concentration of defects c = 0.2, averaged over 8778 realisations, with : ta 
boundary conditions) . The dashed line is based on smoothed data obtained by applying 
a Savitzky-Golay- filter [20] . 

Fig. 5: Effective dynamic exponent %s as a function of time (measured in units of MCA) 
for various atomistic mechanisms and return probabilities. From top to bottom: 
evaporation-condensation (solid), Pda(0 ^ (solid), oc (solid), oc (dotted), 
= 1/2 for Z = ±1 (dashed). Data was smoothed using a Savitzky-Golay-filter. Steps 
with 256 sites, at ksT/e — 1 were simulated. 

Fig. 6: /fc vs. k/M (in units of 27r) for different transport mechanisms: evaporation- 
condensation (plus), terrace diffusion with Pda(0 °^ (cross) and step diffusion with 
Pda(0 = 1/2 for i = ±1 (asterik). The lines correspond to the power-law behaviour ex- 
pected from Langevin theory: Ik oc k'^ (sohd), Ik oc k"^ (dashed) and Ik oc k^ (dotted), 
respectively. Steps with 64 sites, at ksT/e = 0.8, were simulated. 
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